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ABSTRACT 


The  difficulties  encountered  when  implement at ing  ortho- 
graphic reprojection  of  3-dimensional  image  data  onto  a 
2-dimensional  screen  are  considerable.  They  arise  princi- 
pally because  of  the  size  of  data  being  manipulated  and  the 
tendency  for  underlying  or  overlying  structures  to  obscure  a 
clear  view  of  the  desired  image.  In  this  work,  implementa- 
tion was  performed  with  an  orthographic  reprojection  tech- 
nique and  many  heuristic  approaches  were  used  to  resolve 
some  of  the  related  problems. 

Several  coordinate  rotating  algorithms  were  tested  in 
this  work.  Among  them  the  Precalculation  and  Indexing 
method  proved  to  be  the  most  efficient  algorithm. 

Due  to  the  disparity  of  the  viewing- coordinate  grids  and 
the  voxels  of  the  volume  data  after  rotation,  3-dimensional 
interpolation  is  required  for  appling  the  reprojection  tech- 
nique. Several  methods  implementing  linear  interpolation 
have  been  tried.  Interpolation  with  a  cone-shape  kernel  is 
the  most  appropriate  method  in  the  2D  situations  and  can  be 
easily  extended  to  a  sphere-shape  kernel  in  3D  situations. 

The  orthographic  reprojection  method  includes  a  single 
plane  dissection  capability.  Results  on  an  artificial  test 
data  are  collected  using  the  above  algorithms. 

Several  resulting  images  are  included  as  well  as  the 
associated  PASCAL  programs. 
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I.  INTRODUCTION 

A.   GENERAL 

The  representation  of  3 -dimensional ( 3D)  objects  on  a 
2-dimensional (2D)  screen  presents  an  interesting  problem  in 
computer  science  and  computer  engineering.  The  application 
of  this  technology  is  common  in  areas  such  as  computer 
graphics,  CAD/CAE,  and  medical  imaging.  The  study  here  will 
be  concerned  with  a  more  specialized  area  of  showing  3D 
objects  on  2-dimensional  (2D)  screen.  A  variety  of 
techniques  have  been  developed  for  this  purpose. 

A  discrete  set  of  3D  volume  data  consists  of  a  collec- 
tion of  rectangular  parallelepipeds.  These  are  referred  to 
as  volume  elements,  or  more  commonly  voxels.  Each  voxel  has 
an  associated  value  which  represents  the  average  intensity 
of  a  parallelepiped  referred  to  as  the  density  value.  The 
density  value  is  an  integer  from  a  finite  set  (D-j^  .  .  .  D-^^  )  . 
If  D-jnj-n  =0  and  D-mdX,  ~  1  >  we  saY  that  the  image  is  a  binary 
image . 

It  is  useful  to  classify  3D  display  techniques  into 
different  categories.  [Ref.  1] .  The  main  criterion  used  to 
classify  them  is  whether  the  technique  displays  a  whole 
scene  or  a  only  cross- sect ion  of  a  selected  plane  in  the  3D 
volume  data. 

WHOLE  SCENE   DISPLAY  TECHNIQUE The  entire  scene   of  3D 

data  can  be  displayed  by  the  use  of  a  vibrating  mirror.  A 
single  slice  of  a  3D  object  can  be  displayed  onto  a  2D 
screen  at  a  short  time.  If  we  then  display  many  such  slices 
in  a  rapid  succession  and  view  the  screen  through  a 
vibrating  mirror  of  variable  focal  length,  we  can  make  the 
slices  appear  to  be  in  their  correct  spatial  location 
relative  to  each  other. 
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Another  way  of  displaying  the  whole  scene  is  to  project 
it  onto  the  2D  screen.  The  work  studied  here  is  related  to 
this  approach. 

SURFACE   DISPLAY  TECHNIQUE Using   the  surface   display 

technique,  we  are  interesting  in  dipicting  the  appearance  of 
surfaces  of  an  object  in  the  scene.  This  involves 
segmenting  the  original  volume  data  into  an  object  and  a 
background,  for  example,  by  creating  a  binary  image  by  which 
a  1  is  assigned  to  a  voxel  in  the  object  and  a  0  to  a  voxel 
in  the  background.  Various  techniques  such  as  ID  contour 
based,  2D  surface  based,  and  3D  parallelpiped  based  or 
sphere  based  display  techniques  are  examples  of  surface 
displays . 

REPROJECTION  TECHNIQUE There  are  two  kind  of  tech- 
niques in  the  reprojection  method:  Orthographic  reprojection 
and  Radial  reprojection.  [Ref.  2].  Orthographic  reprojec- 
tion is  performed  numerically  in  the  computer  by  summing  the 
values  of  voxels  along  parallel  paths  through  the 
reconstructed  computed  tomography (CT )  scan  image. 

Radial  reprojection  is  an  alternate  scheme  in  which  the 
voxels  of  the  CT  scan  reconstructed  .  volume  are  projected 
onto  a  cylindrical  surface  surrounding  the  object,  rather 
than  orthographically  onto  a  plane. 

The  reprojection  of  volume  image  data  onto  a  2D  screen 
sometimes  cannot  provide  sufficient  information  because  the 
size  of  any  given  structures  is  not  dependent  on  the 
distance  from  the  observer.  Perspective  does  not  exist  in 
the  orthographic  reprojection.  Therefore  it  is  necessary  to 
rotate  the  volume  image  data  to  the  desired  viewing- angle 
before  reprojection.  In  addition  to  the  rotation,  numerical 
dissection  and  numerical  dissolution  techniques  are  intro- 
duced to  show  selected  portions  of  the  reconstructed  volume. 
[Ref.  2]. 
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THE   SCOPE  OF   THE  STUDY The   study   of  techniques   of 

orthographic  reprojection  of  3D  objects  after  rotation  to  a 
selected  viewing- angle  is  the  focus  in  this  work.  The 
computer  technique  for  single  plane  dissection  is  also 
incorporated  here.  Generally  3D  volume  data  consists  of  a 
huge  amount  of  voxels.  Therefore,  an  efficient  algorithm  is 
essential  for  the  manipulation  of  volume  data. 

After  rotation  to  the  selected  viewing-angle ,  each  voxel 
of  the  volume  data  does  not  lie  on  the  grid  of  the  viewing- 
coordinate  system.  In  other  words,  each  voxel  looks  like  a 
polygon  rather  than  a  parallelepiped  along  the  viewing-angle 
from  the  display  screen.  This  will  yield  an  obscuring  effect 
between  voxels  during  projection.  Therefore,  interpolation 
is  necessary  to  resolve  the  contribution  of  voxels  to  each 
pixel  in  the  viewing-plane . 

The  main  emphasis  in  this  work  is  on  the  study  of 
efficient  rotation  algorithms  and  the  interpolation  method. 

In  Chapter  II,  various  rotation  algorithms  are 
discussed,  as  well  as  the  results  obtained  through  the  use 
of  artifical  test  data. 

In  Chapter  III  different  interpolation  methods  using 
various  kernel  functions  are  presented. 

B.   TECHNIQUES  OF  ORTHOGRAPHIC  REPROJECTION 

1 .   Reproj  ect ion  Image  Generation 

The  process  of  reprojection  is  illustrated  in 
Figurel.l.  A  linear  array  of  squares  on  the  viewing-plane 
corresponds  to  the  picture  element s (pixels )  at  the  level  k 
of  this  cross  -  sect  ion .  The  intensity  value  of  the  pixels  in 
the  array  is  the  sum  of  those  voxels  along  the  projection 
path  which  are  perpendicular  to  the  projection  plane.  When 
all  volume  elements  at  the  level  k  have  been  projected  then 
the  voxels  at  level  k+ 1  are  projected. 
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Rotation  of  Volume 


Viewing-Plane 


Reconstructed 
Voxels  at 
Level  K 


Reprojected  Pixels 
at  Level  K 


Figure  1.1    Reprojection  Process. 

The  resulting  2D  array  of  pixel  values  can  be 
rescaled  and  displayed  onto  the  specific  display  device. 
Subsequently  it  is  viewed  as  an  image  on  a  TV  monitor. 

2 .   Numerical  Dissection 

The  dissection  technique  is  used  to  reduce  the 
effect  of  obscuration  and  to  clearly  show  an  internal  struc- 
ture of  the  volume  data.  In  this  technique,  a  plane  of 
dissection  is  highlighted  by  selectively  dimming  the  voxels 
in  front  of  the  plane  and  ignoring  all  the  voxels  behind  the 
plane.  Dimming  is  accomplished  by  replacing  the  value  of 
each  voxel  with  the  product  of  an  original  voxel  density 
value  and  a  constant  less  than  unity.  For  example,  a 
constant  0.1  would  result  in  a  intensity  reduction  to  10  % 
of  its  original  value. 

The  result  is  that  the  internal  structures  at  the 
plane  of  dissection  are  more  clearly  visible.  It  is  also 
easy  to  see  the  spatial  relationship  of  the  object  plane  to 
those  structures  in  front  of  the  plane. 
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3 .   Numerical  Dissolution 

Selective  numerical  dissolution  is  a  process  whereby 
the  relative  contributions  of  selected  voxels  can  be 
decreased  before  rotation.  In  this  manner  obscuring  struc- 
tures are  only  partially  removed  and  this  results  in  a  'see 
through'  effect. 

The  criterion  for  numerical  dissolution  is  based  on 
the  density  difference  between  the  structures  within  the 
reconstructed  volume.  The  density  threshold  is  identified 
empirically  by  first  choosing  an  intitial  threshold  and  then 
adjusting  the  value  to  achieve  the  desired  dissolution 
effect.  Once  the  desired  threshold  is  identified,  numerical 
dissolution  is  accomplished  by  dimming  all  voxels  with 
values  below  threshold  during  reproj ect ion .  The  dimming 
process  is  the  same  as  that  in  the  dissection  process. 
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II.  COORDINATE  TRANSFORMATION 

A.   GENERAL 

1 .   Introduction 

Several  methods  have  been  proposed  for  the  display 
of  3 -dimensional ( 3D)  information  contained  within  a  series 
of  parallel  Computed  Tomography (CT )  cross  sections.  This  is 
referred  to  as  a  3D  or  volume  reconstruction  problem.  The 
reconstructed  volume  data  will  be  displayed  on  a 
2-dimensional ( 2D)  TV  screen  using  reproj ect ion ,  dissection, 
and  dissolution  techniques.  Although  the  problem  of  obscu- 
ration exists  in  the  reprojected  image,  it  is  not  as  severe 
as  in  the  case  of  conventional  radiographs  because  the 
reconstructed  volume  can  be  viewed  from  different  viewing- 
angles .  Before  reprojection  the  object  can  be  mathematically 
oriented  to  a  desired  viewing- angle .  The  choice  of  desired 
viewing-angles  will  enable  radiologists  to  discover  the 
important  information  from  the  image. 

The  objective  of  this  chapter  is  to  present  an 
example  of  3D  coordinate  transformation.  Because  digital 
image  data  is  comprised  of  a  huge  amount  of  3D  data  points, 
our  rotation  algorithm  must  be  efficient.  Also,  coordinate 
rotation  invokes  time-consuming  floating-point  operations  in 
the  computers.  An  inefficient  algorithm  will  prove  unwieldy 
and  lead  to  an  undesirably  long  processing  delay  before 
results  are  available.  In  the  final  part  of  this  chapter, 
hardware  alternatives  to  the  software  design  which  will 
implement  the  rotation  algorithm  developed  in  this  chapter 
will  be  considered. 
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Center  of  Rotation 


The  center  of  rotation  is  arbitrary,  although  it  is 
generally  chosen  to  be  near  the  center  of  the  object  (Figure 
2.1).  This  is  done  so  that  after  rotation  the  reconstructed 
image  will  not  be  shifted  too  much  out  of  the  effective 
viewing-window  defined  on  the  screen. 

Z 


Center  oil  Volume 
Rotation  about 
X-Axis 


Projected 
of  an  Object 


Potation  about 
Y-Axis 


Ce  iter  of  Scr^  en 


X 


Full  Screen  Size 


(A)  (B) 

Figure  2.1    Center  of  Rotation(A)  and   Viewing-Plane (B ) . 

3 .   Rotation  Types 

Three  kinds  of  rotation  are  used  for  this  work: 
Rotation  about  X-axis  (Elevat ional ) ,  rotation  about  Y-axis 
(Horizontal),  and  rotation  about  both-axes(a  combination  of 
rotation  about  the  X-axis  followed  by  rotation  about  the 
Y-axis).  These  three  kinds  of  rotation  will  be  enough  to 
provide  a  desired  viewing-angle  from  any  directions.  Figure 
2.2  illustrates  examples  of  elevational  and  horizontal  rota- 
tions. Figure  2.2(A)  shows  +90  degrees  rotation  about  X-axis 
and  Figure  2.2(B)   shows  +90   degrees  rotation  about  Y-axis. 
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(A)90   Rotation    (B)Original  Image    (C)90   Rotation 
about  X-Axis  about  Y-Axis 


Figure  2.2 


Direction  of  Rotation 


Rotation  about  both-axes  is  performed   by  the  combination  of 
elevational  and  horizontal  rotations  sequentially. 

4 ,   Coordinate  System 

Figure  2.3  shows   the  coordinate  system  we   will  use 
throughout  this  work. 

A  coordinate  system  is  selected  according  to  the  right 
hand  rule.  As  shown  in  Figure  2.3,  the  horizontal  line 
represents  the  X-axis,  the  vertical  line  represents  the 
Y-axis,  and  depth  is  labelled  with  the  Z-axis. 
Obj ect -  coordinates  are  those  we  define  relative  to  the 
object  during  image  sampling,  and  viewing- coordinates  are 
those  defined  for  displaying  the  image.  The  viewing- 
coordinate  system  always  corresponds  to  the  coordinate  of  a 
display  screen.  Thus  the  orientation  of  the  object  with 
respect  to  the  viewing- coordinates  defines  the  rotation 
angle . 
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♦»  X 


X,Y,Z  :Obiect-Coordinate 
X» ,¥• ,Z« : Viewing-Coordinate 


Figure  2.3    Coordinate  System 


B.   COORDINATE  TRANSFORMATION 


1 .   Matrix  Multiplication 


Usually  a  coordinate  transformation  is  done  by 
multiplying  the  original  coordinate  by  a  transformation 
matrix.  If  we  use  matrix  notation  this  will  result  in  a 
vector  transformation  from  an  R3  space  to  another  R3  space. 
We  can  define  a  function  T:  R3  -->  R3  by  T(X)  =  AX  where  A 
is  a  transformation  matrix  and  X  is  a  input  vector.  Observe 
that  if  X  is  a  3x1  matrix  and  A  is  3x3  matrix,  then  the 
product  AX  is  also  a  3x1  matrix.  Thus  T  maps  R3  into  R3 . 
This  kind  of  linear  transformation  is  called  'matrix  trans- 
formation'. Algebraically  a  3D  to  3D  transformation  matrix 
requires  3  basis  functions.  Suppose  we  want  to  rotate  a  3D 
vector  X  =  xx  +  y£  +  z£.  The  transformation  matrix  for 
•this  operation  should  be  a  3x3  matrix.  In  matrix  notation  we 
write 
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'ax  +  by  +  cz 
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^gx  +  hy  +  iz, 
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If  the  determinant  of  A  is  1,  it  produces  a  pure  rotation 
about  the  origin.  Before  considering  rotation  about  an  arbi- 
trary axis,  two  special  cases  are  examined:  rotation  about 
the  X-axis  ( elevational )  and  rotation  about  the 
Y- axis (horizontal ) . 

2 .   Rotation  about  X-Axis 

In  this  case  the  rotation  matrix  will  have  zero  in 
the  first  row  and  first  column  except  a  unit  value  on  the 
main  diagonal.  The  other  terms  are  rotation  dependent  and 
will  be  determined  by  considering  the  specific  rotation 
angle  of  a  point.  Figure  2.4  shows  us  this  situation  geome- 
trically. Rotation  about  only  one  axis  is  the  same  as  2D 
transformation  because  the  rotational  axis  is  fixed  in 
space.  Here  the  rotation  angle  is  defined  to  be  positive  in 
the  counter-clockwise  direction  when  a  vector  rotates. 

.a. 
*  In  this  figure,  P  is  an  unit  vector. 


f(x\y\z) 


~  -|-^»ip(x,y,z) 


Figure  2.4 


Rotation  about  X-Axis 


In  the  Fig  2.4,  the  X,  Y,  and  Z  represent  the  object- 
coordinates  and  X# ,  Y» ,  and  Z»  the  viewing- coordinate  system 
respectively.  Thus  our  matrix  transformation  fomula  is 
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P(x,y,z)  =  =  =  >   P»(x»,y,z«) 


x» 
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cos  8 
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sinO 
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f       X 
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sing 

y 

cos  6 
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In  figure  2.4,   we  let  x  =    x, 
rotation  we  have 


x       x 

ycosQ  -  zsin6 
ysinQ  +  zcos  Q 


y  =  cos<p, 


z  =  si 


ru^>. 


After 


x#  =  x  , 

y •  =  cos(  B   +  d>)    =  cos9cos0-  sinQsin^'  =  ycos  #  -  zsin£, 

z»  =  sin(^  +  d>)    =  sinQcos^  +  cosgsin^  =  ysin£  +  zcosg. 

This  result  is  just  the  matrix  multiplication,  the 
product  of  a  transformation  matrix  and  a  vector  in  the 
ob j  ect -  coordinate . 

3 .   Rotation  about  Y-Axis 

Rotation  about  the   y-axis  can  be  shown   in  the  same 
manner  as  rotation  about  the  X-axis  above.  (See  Figure  2.5) 


X 


y 

z* 


P(x,y,z)  =  =  =  >  P»(x»,y,z») 

_x 


cos  9    0    sinB 

0     1       0 

-sin£)     0    cosB 
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<zs 


xcos  6  +  zsin  0 

y 

- xsinQ  +  zcos  Q 


S 


4 .   Rotation  about  Both- Axis 

Rotation  about  an  arbitrary  axis  requires  a  more 
complex  operation.  But  we  can  restrict  ourselves  to  two 
rotation  directions  for  adequate  display  purposes.  This 
combination  of  two  directional  rotations  can  be  considered 
as  a  cascaded  operation,  first  about  the  X-axis  and  then 
about  the  Y-axis  or  vice  versa.  If  rotation  about  the 
X-axis,  then  rotation  about  Y-axis  is  chosen,  then  the 
rotation  matrix  will  be 
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z  ^  z 


Z  4 


i_.2p(x>y>z) 


Figure  2.5    Rotation  about  Y-Axis 
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This  matrix  could  not  be  applied  in  the  reverse  direc- 
tional rotation,  because  a  matrix  multiplication  AB  is 
generally  not  equal  to  BA.  Therefore  we  have  to  be  careful 
in  our  selection  of  matrix  multiplication  sequences. 

C.   PROGRAMMING  METHODS 

General  procedures  for  the  image  display  are  as  follows: 
(1)  Set  up  the  obj ect -  coordinate  value,  (2)  Calculate  the 
corresponding  viewing- coordinate  value  by  coordinate  trans- 
formation, (3)  Access  a  pixel  and  place  it  in  the  corre- 
sponding viewing- coordinate ,  (4)  Perform  projection  and 
display . 

The  following  algorithms  use  different  methods  for  coor- 
dinate transformation.  Because  coordinate  transformation  is 
the  central  issue  in  this  problem  the  program  efficiency 
will  mainly  depend  on  it. 
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1 .   Direct  Matrix  Multiplication 

Suppose  we  want  to  rotate  a  3D  object.  Every  data 
point  in  the  obj ect -  coordinate  has  x,  y,  z  coordinates  which 
form  a  3D  position  vector.  To  rotate  this  object  we  have  to 
multiply  the  position  vector  by  a  3x3  transformation  matrix. 
As  shown  in  Table  1,  the  transform  of  each  vector  requires  9 
real  multiplications  and  6  real  additions.  If  one  dimension 
of  the  volume  is  N,  the  total  calculation  reaches  9N3  multi- 
plications and  6N3  additions.  The  flow  chart  for  imple- 
menting this  method  is  shown  in  Figure  2.6.  In  a  PASCAL 
environment,  one  record  consists  one  plane  and  each  record 
is  accessed  sequentially.  Whenever  an  image  plane  is 
accessed  the  program  rotates  the  coordinates  of  pixels 
contained  in  the  plane  one  by  one.  Therefore  the  subroutine 
'Rotation'  is  called  N3  times. 


Define  Viewing-Plane 


Open  Data  File 


Get  the  Direction  of 
Rotation  and  the  Angle 


Calculate  Elements  of 
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Figure  2.6    Flowchart  of  the  Direct  Matrix  Multiplication 
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2.   Decomposition  with  Fixed  Coordinate 

Let  us  consider  a  rotation  about  only  one  axis.    As 

shown  in  Equation   2.1  the  transformation  matrix  includes  4 

zeros  and   a  unit  value  in   the  main  diagonal  so  the  center 
axis  of  the  rotation  will  not  be  changed 


x« 
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1      0 
0   cos  0 

[0      sin  0 


0Vxs 


sin£ 
cos  QJ 


y 
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ycos  9  -  zsin  0 
zcos  9   +  zcos  9 


X'l  +  y  •  0  +  z  •  0 

x*0  +  y.cos  ^  -  z«sin9 

x»0  +  y»sin9  +  z*cos©^ 

(2.2) 


(2.1) 


Equation  2.1  shows  an  example  of  rotation  about  the 
X-axis  using  the  direct  matrix  multiplication  method. 
Equation  2.2  suggests  a  more  efficient  implementation 
method.  If  the  matrix  multiplication  is  simplified  to  3 
equations  as  shown  in  Equation  2.2,  it  results  in  an  economy 
of  4  multiplications  and  4  additions  for  each  pixel 
transformation . 

Another  point  in  comparing  the  two  equations  is  that 
the  first  row  of  the  transformation  matrix  consists  of  only 
a  unit  value  and  two  zeros.  This  row  maps  a  incoming  vector 
into  itself.  (This  is  the  fixed  coordinate.)  The  other  two 
coordinate  values  will  be  changed  according  to  Equation  2.2. 
This  method  is  represented  in  the  flow  chart  shown  in  Figure 
2.7.  The  inner  loop  does  not  include  a  call  for  a  matrix 
multiplication  operation,  but  the  outer  two  loops  call  the 
subroutine  'Rotation'  which  includes  matrix  multiplication 
operations.  Therefore  the  total  number  of  calculations 
required  are  4N2  multiplications  and  2N2  additions. 

This  method  of  decomposition  with  fixed  coordinates 
is  not  very  helpful  for  the  case  of  arbitrary  directional 
rotation.    The   transformation   matrix    of   the   arbitrary 
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directional  rotation  has   only  one  zero  in  the   first  row  so 
the  economy  of  calculations  will  be  very  small. 
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Figure  2.7    Flow-Chart  of  Decomposition  with  Fixed 
Coordinate (Rotation  about  One-axis). 


3 .   Precalculation  and  Indexing  Method 

Another  efficient  method  makes  use  of  indexing 
instead  of  multiplication.  This  method  is  represented  in  the 
following  equation 


'x-s 


y 


z  , 


ell  cl2   cl3 
c21  c22   c23 


n  r   s 

X 


c31  c32   c33 


\ 


/ 


kz/ 


cllx  +  cl2y  +  cl3z 
c21x  +  c22y  +  c23z 
c31x  +  c32y  +  c33z 


(2.3) 
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The  elements(C'j  )  of  the  transformation  matrix  are  fixed 
after  the  rotation  angle  is  specified,  and  the  x,  y ,  and  z 
values  of  the  obj ect -  coordinate  iterate  individually.  We 
can  calculate  the  multiplication  of  each  coordinate  and 
corresponding  matrix  element  (  CjjX|<)  before  accessing  pixels 
of  volume  and  then  store  those  values  in  an  array.  After 
this ,  by  indexing  the  array  according  to  the  incoming 
obj ect- coordinate  and  adding  these  components  together,  we 
can  calculate  the  viewing- coordinate  value.  Implementing 
this  idea  can  allow  us  to  avoid  the  time-consuming  floating- 
point multiplications.  It  is  represented  in  the  flowchart 
of  Figure  2.8.  As  shown  in  the  figure  the  Indexing  method 
does  not  use  a  'Rotation'  subroutine. 

The  array  contents  also  can  be  calculated  by  succes- 
sive additions  as  shown  in  the  program  ' Rotat ion_3_Dim ' . 
When  the  multiplication  of  a  coordinate  and  matrix-element 
is  calculated,  a  coordinate  value  always  increases  by  one 
unit.    For  example,    X^+i:  =  X^  +  1,    so  the  multiplication 

(CjjXK)  becomes  C^X^  CUXK+  ^U  '  CliXK^=  CHXK+1  +  C„  •••etc. 
Only  the  first  term  requires  multiplication  of  4  pairs  of 
elements.  The  total  number  of  additions  required  is  4N 
+  2N2.  The  4N  additions  are  the  effective  multiplications  of 
the  objective  coordinate  and  transformation  matrix  element. 
The  remaining  2N2  additions  result  from  calculations  of  the 
viewing- coordinate  value  by  means  of  double  iteration  loops. 
Because  all  floating-point  multiplications  are 
avoided  by  indexing  and  adding,  this  method  results  in  large 
processing  time  savings.  The  total  number  of  calculations 
required  is  shown  in  Table  1.  Compared  to  the  other  two 
methods,  Table  1  indicates  that  this  method  is  the  most 
efficient  of  the  three.  But  this  method  requires  some 
working  space  in  memory  for  storing  the  precalculated 
values.  In  a  VAX-11/750  floating-point  notation  (where  a 
floating   point   value   occupies  4   bytes)    12N   bytes   are 
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required  for  a  rotation  about  one  axis,  and  24N  bytes  are 
required  for  rotation  about  both  axes.  This  is  a  reasonable 
memory  requirement  considering  the  processing  time  savings 
realized  by  using  this  method. 
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Figure  '2.8    Flow  Chart  of  Precalculation  and  Indexing 
Method (Rot at  ion  about  one-axis). 


D.   PROGRAM  IMPLEMENTATION 

The  program  'Rot_3_Dim'  which  implements  the 
Precalculation  and  Indexing  method  is  included  in  Appendix 
A.  In  this  PASCAL  program,  'Objected'  and  'Viewed'  repre- 
sent the  obj ect -  coordinate  and  the  viewing- coordinate 
respectively.    Data     used   to   test   the   program    is  an 
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artificially  created  double  pyramid(  N  =  64  ).  This  3D  data 
is  shown  in  Figure  2.9  which  portrays  the  shape  of  a  peak 
connected   double  pyramid.    The  density   of   the  test   data 


Figure  2.9    Double  Pyramid. 
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Viewing-...., 
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Full  COMTAL  Screen 


(A)  (B) 

Figure  2.10    Test  Data(A)  and  Display  Buffer(B). 

increases  toward  the  center  of  pyramid.  The  object  coordi- 
nate is  defined  so  that  the  center  of  object  is  located  at 
the  center  of  the  coordinate.    The  size  along  one  dimension 
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of  the  data  is  64,  but  the  program  needs  111x111  elements 
size  of  a  viewing-plane .  In  the  worst  case  one  side  of  the 
reconstructed  image  can  have  the  diagonal  size  of  the 
original  image (Figure  2.10). 

E.   NUMBER  OF  CALCULATIONS  AND  PROCESSING  TIME 

The  program  'Rot_3_Dim'  which  manipulates  our  test  data 
is  run  on  the  DEC  VAX- 11/750  computer.  Results  are 
displayed  on  the  COMTAL  image  processing  system  which  itself 
has  a  PDP-11/23  processor.  Table  1  shows  a  comparison  of 
program  efficiency  and  execution  time  among  these  3  methods. 


TABLE  I 
Table  1.  Comparision  of  Methods 


\.  Direction 
Methods   \. 

Rotation  about 
One-Axis 

Rotation  about 
Both- Axis 

NO  of  Cal. 

Time  Taken 

NO  of  Cal. 

Time  Taken 

Direct  Matrix 
Multiplication 

9N3  Mults 
6N3  Adds 

over  25 
sees 

9N3  Mults 
6N3  Adds 

over  10 
mins 

Decomposition 
with  Fixed 
Coordinate 

4N2  Mults 
2N2  Adds 

25  Sees 

4N2+4N3 
Mults 

2N2+2N3 
Adds 

10  Mins 

Precalculation 
and  Indexing 

4  Mults 
4N+2N2  Adds 
(16N  Bytes 
Storage ) 

15  Sees 

8  Mults 
9N3  Adds 
(32N  Bytes 
Storages ) 

25  Sees 

-  N  :  The  size  of  dimension. 
All  real  value  calculation. 

F.   HARDWARE  REALIZATION 

Modern  computers  have   only  one  CPU.   In   this  case  each 
line  of   code  must   be  executed   sequentially.   But   digital 
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image  data  is  composed  of  a  huge  number  of  pixels.  Therefore 
a  special  purpose  computer  having  parallel  processing  capa- 
bility is  becoming  popular.  Here  two  simple  hardware  designs 
are  suggested  to  realize  the  above  algorithms  of  coordinate 
rotation.  This  hardware  decreases  the  burden  of  iterative 
matrix  multiplication  in  the  coordinate  rotation  and 
releases  the  CPU  to  do  more  important  task.  Also  it  will 
increase  the  speed  of  the  image  rotation.  These  two  design 
schemes  use  different  components.  Figure  2.11  illustrates 
the  first  hardware  design  scheme.  The  first  design  consists 
of  3  adders  and  9  floating-point  multipliers.  An  incoming  X 
coordinate  is  multiplied  with  the  first  column  of  the  trans- 
formation matrix,  a  Y  coordinate  the  second  column  , and  a  Z 
coordinate  the  third  column  respectively.  Then  the  multi- 
plied coordinates  and  matrix  elements  are  added  to  make 
viewing- coordinate  values  as  shown  in  Figure  2.11. 

Figure  2.12  shows  the  second  hardware  design.  The  second 
one  consists  of  3  adders  and  9  /  4N  bytes  memories  instead 
of  multipliers.  Each  set  of  the  memory  has  the  capacity  of 
4N  bytes  and  individual  memory  element  has  4  bytes  width. 
The  precalculated  coordinates  and  matrix  elements  will  be 
stored  in  the  memories.  These  components  are  indexed 
according  to  the  incoming  obj ect -  coordinate  values,  and  then 
added  to  calculate  viewing- coordinate  values.  It  is  clear 
that  an  indexing  is  faster  than  a  floating  point  multiplica- 
tion. As  may  be  expected  there  is  a  trade-off  between  cost 
and  speed. 


30 


X 


Yo— 


=<D*Z- 


Figure  2.11    First  scheme(  9  /  Floating-Point  Multipliers) 
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Figure  2.12    Second  Scheme(  9  /  4N  Bytes  Memories  ). 
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III.  INTERPOLATION 

A.  INTERPOLATION  IN  IMAGE  PROCESSING 

Interpolation  is  the  process  of  estimating  the  interme- 
diate values  of  a  continuous  signal  from  discrete  samples. 
Interpolating  is  used  extensively  in  digital  image 
processing  to  magnify  or  reduce  images  and  to  correct  a 
spatial  distortion.  In  3 -dimensional ( 3D)  digital  image 
processing,  each  volume  element (voxel )  value  represents  the 
intensity  value  of  a  rectangular  parallelepiped.  This 
rectangular  parallelepiped  is  not  necessarily  a  cube,  but  it 
is  assumed  to  be  a  cube  that  has  equal  length  edges  for  our 
purposes  here.  A  digital  intensity  value  is  generated  by 
analog- to-digital  conversion  of  a  continuous  intensity 
signal.  Therefore,  a  digital  data  produced  in  this  way  can 
involve  a  certain  amount  of  aliasing  or  blurring  in  the 
digitization  process.  But  this  effect  is  so  small  that  we 
can  ignore  it.  An  exact  image  restoration  will  thus  be 
possible  if  an  appropriate  interpolation  method  is  applied. 
The  goal  of  the  interpolation  of  an  encoded  medical  image  is 
to  reconstruct  the  original  image  as  closely  as  possible. 
Because  of  the  great  amount  of  data  associated  with  a 
digital  image,  an  efficient  interpolation  algorithm  is 
essential . 

B.  INTERPOLATION  METHOD 

1 .   One -Dimensional  Interpolation 

An  interpolation  kernel  function  is  a  special  type 
of  approximating  function.  A  fundamental  property  of  an 
interpolation   kernel  function   is   that   it  must   have   the 
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S(X)  =  SQUARE  ******  square,   *  :  convolution 
n  times 
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Figure  3.1    Interpolation  Functions ( ID) 
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sampled  data  values  at  the  interpolating  knots  or  sampled 
points.  In  other  words,  if  F  is  a  sampled  function  and  G  is 
the  corresponding  interpolating  function,  then  G(Xk)  -  F(Xk) 
whenever  Xk  is  an  interpolation  knot.  For  equally  spaced 
data  samples,  many  interpolation  functions  can  be  written  in 
the  form 

G(x)  =£CRU(-^)  (3.1) 

where  h  represents  the  sampling  increment,  the  Xj^'s  are  the 
interpolation  knots, and  U  is  an  interpolation  kernel.  The 
interpolation  kernel  in  Equation  3.1  converts  discrete 
data  into  a  continuous  function  by  an  operation  similiar  to 
the  convolution.  Mathematically  developed  interpolation 
kernel  functions  include  nearest  neighboring,  linear,  cubic 
spline,  cubic  convolution,  piecewise  linear,  and  sine  func- 
tion. A  one-dimensional ( ID)  version  of  these  interpolation 
functions  are  illustrated  in  Figure  3.1.  [Ref.  3].  The 
first  function  in  Figure  3 . 1  is  a  sample  and  hold,  nearest 
neighboring,  or  replication  interpolation.  The  second  func- 
tion is  a  square  convolved  with  another  square  which  results 
in  a  triangle  or  linear  interpolation.  The  third  function 
is  the  convolution  of  three  squares  or  a  cubic  spline.  The 
fourth  function  is  an  arbitrary  piecewise  linear  function 
which  illuminates  the  arbitrariness  of  the  interpolation 
kernel.  The  fifth  function  is  a  cubic  convolution  kernel 
which  is  composed  of  piecewise  cubic  polynomials  defined  on 
the  subintervals  (-2,-1),  (-1,0),  (0,1),  and  (1,2). 
[Ref.  3].  The  last  one  is  the  sine  function  which  provides 
an  exact  reconstruction.  The  use  of  the  sine  function 
(sinX/X)  as  an  interpolant,  which  has  negative  valued  side 
lobes,  requires  an  infinite  number  of  terms.  Realization  is 
difficult.  difficult  in  the  realization.  Polynomials  of 
order  two  or  greater  can  also  be  employed  as  an  interpola- 
tion kernel  function.   The  amplitude  spectra  of  the  nearest- 
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Figure  3.2    Amplitude  Spectra  of  Interpolation  Functions 


(A)Nearest  Neighboring 


(B)Linear  Interpolation 


Reconstructed  signal 


(C)Sinc  Function 

Figure  3.3    One-Dimensional  Interpolation  Process. 

neighboring,  linear  interpolation,  and  cubic  convolution 
interpolation  kernels  are  shown  in  Figure  3.2  for  frequency 
from   0   to   4TT/h   (where    h   is   the   sampling   interval). 
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[Ref.  3].  The  response  of  an  ideal  interpolation  kernel  is 
a  unit  step  function  in  the  frequency  domain.  In  image 
data,  the  loss  of  high  frequency  information  causes  the 
image  to  appear  blurred.  On  the  other  hand,  deviation  from 
the  ideal  spectrum  beyond  the  shaded  area  contributes  to 
aliasing.  One-dimensional  interpolation  examples  with 
several  interpolation  kernels  are  performed  in  a  fashions 
shown  in  Figure  3.3.  [Ref.  4],  Interpolation  kernels  have 
a  significant  impact  on  the  numerical  behavior  of  the  inter- 
polated functions.  Because  of  their  influence  on  accuracy 
and  efficiency,  it  is  necessary  to  select  carefully  an 
appropriate  interpolating  kernel  for  an  image  processing. 

2 .   Two -Dimensional  Interpolation 
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(A)Piecewise  Linear  Interpolation  (B)Bilinear  Interpolation 
Figure  3.4    Two-Dimensional  Linear  Interpolation. 

A  two-dimensional ( 2D)  interpolation  is  accomplished 
by  two  seperate  ID  interpolations  with  respect  to  each  coor- 
dinate. It  should  be  performed  along  separable  orthogonal 
coordinates  of  the  continuous  signal.  A  2D  linear  interpo- 
lation kernel  function  is  an  example  of  an  orthogonally 
separable  interpolation  function:  G(x,y)  =  G(x)-G(y). 
Figure  3.4   shows  two  examples   of  the   linear  interpolation 
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method  in  2D  situations.  [Ref.  4].  Example  (A)  shown  in 
Figure  3.4  is  performed  in  a  piecewise  fashion.  In  region  I 
of  example  (A),  points  are  linearly  interpolated  in  the 
plane  defined  by  pixels  A,B,C,  while  in  region  II,  the 
interpolation  is  performed  in  the  plane  defined  by  pixels 
B,C,D.  The  continuous  bilinear  interpolation  is  shown  in 
example  (B).  It  is  done  by  linearly  interpolating  points 
along  separable  orthogonal  coordinates  of  the  continuous 
image  signal . 

C.   PROBLEMS  IN  3 -DIMENSIONAL  DATA  PROJECTION 

A  3D  data  is  a  collection  of  discrete  values  resulting 
from  equalspace  sampling.  This  data  consists  of  small  iden- 
tical parallelepipeds  (voxels)  divided  by  three  sets  of 
planes  parallel  to  the  X,  Y,  and  z  axes.  Each  voxel  has  a 
sample  value.  It  is  referred  to  as  the  intensity  value  of  a 
voxel.  For  clarity  of  discussion  the  voxel  value  is  located 
on  the  grid  of  the  obj ect -  coordinate  system.  Figure  3.5(A) 
shows  the  relationship  between  the  voxels  and  the  object- 
coordinate  system  before  rotation.  Orthogonal  projection 
along  the  Z-axis  is  performed  by  integrating  the  intensity 
value  of  voxels  along  a  line  parallel  to  the  Z-axis (  any 
point  P  lying  in  one  of  the  X-Y  plane  in  the  object- 
coordinate  can  not  be  obscured  with  another  point  on  the 
same  plane).  The  obj ect -  coordinate  grids  exactly  correspond 
to  the  viewing- coordinate  grids. 

After  rotation,  the  obj ect -  coordinate  grid  does  not  share 
the  same  location  as  the  viewing- coordinate  grid  as  shown  in 
Figure  3.5(B).  The  coordinate  value  of  a  voxel  on  the  grid 
of  the  obj ect- coordinate  will  be  transformed  to  its  viewing- 
coordinate  values .  These  values  may  not  be  integers  which 
are  the  viewing- coordinate  grids.  This  is  the  reason  why  we 
have  to  interpolate  intermediate  intensity  values  at  the 
viewing- coordinate  grids. 
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»Y,Y 


Viewing- 
Coordinate 


(A) 


(B) 


Figure  3.5    Relationship  between  3 -Dimensional  Data  and 
Viewing- Coordinate  Before  Rotation(A)  and  After  Rotation(B) 


To  reconstruct  the  exact  object,  a  high  order  interpola- 
tion kernel  such  as  sine  functions  will  produce  better 
results.  But  using  a  high  order  kernel  requires  many  calcu- 
lations and  sometimes  it  exceeds  reasonable  computer  capa- 
bility. If  we  can  achieve  an  appropriate  result  with  a  low 
order  kernel  function,  We  can  realize  a  trade-off  between 
fidelity  and  processing  time.  For  that  reason,  only  the 
linear  interpolation  kernel  in  the  2D  and  3D  cases  will  be 
used  in  this  work. 

D.   APPLICATION  OF  LINEAR  INTERPOLATION 

1 .   One -Dimensional  Interpolation 

A  linear  interpolation  function  will  yield  lower 
quality  results  compared  to  that  of  the  high  order  interpo- 
lation function  .  An  intermediate  interpolated  point  is 
calculated  with  the  nearest  two  sample  points.  Linear 
interpolation  includes  an  assumption  that  the  two  adjacent 
sample  points  have  a  linear  relationship. 
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Let  us  consider  a  linear  interpolation  example  of  a 
signal  after  rotation  in  a  ID  case.  Figure  3.6  illuminates 
this    situation. 


&  A     Object -Coordinate 


x±xl*    x2     „•    x3      _.  xU         .  x5 
x2  x3  x4 


8    Rotation 


X " Viewing-Coordinate 


Figure  3 . 6 


Projection  Process  after  Rotation  in 
One -Dimension. 


Suppose  an  arbitrary  function  is  rotated  in  the  counter- 
clockwise direction  by  an  amount  Q  degrees.  In  this  case  the 
sample  points  (PI , P2 , P3 , • • • )  will  be  projected  onto  the 
viewing-plane  (0-A*).  A  projected  viewing- coordinate  value 
at  the  grid  of  viewing-plane  will  have  nonintegers.  In  the 
above  figure,  xl,x2,x3,#,#  are  projected  onto  viewing- 
coordinate  values  with  the  sample  values  (PI ,P2 , P3  ,  • • • ) . 
Therefore  we  have  to  determine  the  intensity  values  at  the 
the  grids  of  obj ect -  coordinate .  Because  the  adjacent  two 
sample  values  have  a  linear  relationship,  we  can  determine 
the  intermediate  values  from  these  two  adjacent  sample 
points  (Figure  3.6) 

PI-  =  PI  *  (x2  -  xl«)  +  P2  *  (xl«  -  xl)  (3.2) 
P2»  =  P2  *  (x3  -  x2«)  +  P3  *  (x2«  -  x2 )  (3.3) 
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As  shown   in  Equation   3.2  and   3.3,   each   interpolation 
requires  2  multiplications,  1  addition,  and  2  subtractions. 

2 .   Two -Dimensional  Interpolation 

The  two-dimensional  interpolation  is  an  extension  of 
a  ID  interpolation  because  the  signal  function  is  defined  in 
the  two  orthogonal  coordinates  as  G(x,y)  =  G(x)G(y).  Here 
three  algorithms  are  tested  with  a  sinusoidal  2D  function 
having  642  elements:  F(x,y)  =  127cos(x+y)  +  127.  For  this 
example  an  identical  sampling  interval  in  both  X,  and  y 
coordinates  are  used.  The  center  of  rotation  is  always 
taken  to  be  the  center  of  the  coordinate  systems.  The  rota- 
tion is  applied  to  the  function  which  is  written  in  object- 
coordinate  (X  ,Y)  .  After  rotation  the  sample  points  of  the 
original  function  no  longer  coincide  with  the  grid  of 
viewing- coordinates .  It  is  necessary  to  interpolate  the 
signal  value  at  the  grids  of  the  viewing- coordinate . 


Center  of 
Rotation 


Object-Coordinate 

X  Viewing-Coordinate 


Figure  3 . 7 


Two-Dimensional  Rotation. 
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a.   Interpolation  by  Distributing   the  Sample  Values 
to  the  Grids  of  the  Viewing-Coordinate 

Since  the  data  values  are  samples  on  the  grids 
of  the  obj ect -  coordinate ,  interpolation  of  the  intermediate 
signal  should  be  done  in  the  obj ect -  coordinate .  But  for  the 
practical  reason  of  data  storage  the  interpolation  by 
density  distribution  in  the  viewing- coordinate  is  attempted 
first.  In  this  case,  a  linear  interpolation  kernel  on  the 
viewing- coordinate  is  assumed  similiar  to  the  kernel  on  the 
obj ect -  coordinate .   The  algorithm  for  this  method  is 

(1)  Consider  one  obj ect -  coordinate  grid. 

(2)  Transform  it  into  a  viewing- coordinate  value  by 
coordinate  rotation. 

(3)  Distribute  a  pixel  value  inversely  proportional  to 
the  distance  from  the  nearest  4  surrounding  viewing- 
coordinate  grids.   (see  Figure  3.8(a)) 


Pixel  Val 


Obj ect -Coordinate 
Assumed  Coordinate 
Parallel  to  Viewing-  Y,  y 


Coordinate   P2 


Viewing -Coordi 


P3(ix+I,iy+1Y) 


►  X 
pr(ix,iy)    \P4(ix+l,iy) 


(A) 


(B) 


Figure  3.8    Interpolation  Kernel  over  Viewing- Coordinate (A) 
Sampling  Value  Dis t ribut ion(B ) . 
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pi-  =  [  Po  *  (l-(vx-ix))] [l-(vy-iy)] 

P2«  =  [  Po  *  (1- (vx-ix) )] [vy-iy] 

P3«  =  [  Po  *  (vx-ix)] [ (vy-iy)] 

P4-  =  [  Po  *  (vx-ix)] [1- (vy-iy)] 

Where  Pl»,P2»  ,P3*  ,P4-  are  distributed  components  of 
the  4  neighbors  of  the  original  sample  Po  located  at 
(vx,vy).  The  displayed  image  thus  obtained  is  shown  in  the 
Figure  3.9.  This  figure  indicates  that  the  error  is  propor- 
tional to  the  rotation  angle  and  reaches  a  maximum  at  45 
degrees  rotation.  As  shown  in  the  ID  interpolation  case, 
the  pixel  value  of  the  rotated  obj ect- coordinate  grid  will 
be  projected  orthogonally  onto  the  viewing-plane  (x*-y# 
plane).  Even  though  the  object  coordinate  grids  are  evenly 
spaced,  the  projected  grids  become  compressed  on  the 
viewing-plane.  But  this  method  assumes  the  evenly  spaced 
grids  regardless  of  a  rotation  angle.  Since  the  amount  of 
error  introduced  depends  on  the  rotation  angle,  it  is 
impractical  to  apply  it  in  the  real  display. 

b.   Interpolation  of   the  Signal  from  the   Nearest  4 
Samples  at  the  Grids  of  the  Obj ect -Coordinate 

In  this  algorithm,  interpolation  is  performed  in 
bilinear  fashion  on  the  obj ect -  coordinate .  The  algorithm 
for  this  method  is 

(1)  Consider  one  viewing- coordinate  grid. 

(2)  Transform  it  into  obj ect -  coordinate  value  by  an 
inverse  coordinate  rotation. 

(3)  Perform  the  interpolation  -from  its  nearest  4 
surrounding  samples  by  the  bilinear  method. 

This  is  a  correct  linear  interpolation.  However  problems 
with  this  algorithm  are  (1)  We  have  to  rotate  all  viewing- 
coordinate  grids.    (2)    All  input  data  should   be  accessed 
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(A)  Original  Image 


(B)  15  Degrees  Rotation 


(C)  30  Degrees  Rotation       (D)  45  Degrees  Rotation 


Figure  3.9    Displayed  linage  with  the  Interpolation 
by  Distributing  the  Sample  Values  to  the  Grids 
of  the  Viewing- Coordinate . 


simultaneously   in  computer   memory,    which  requires   large 
working  space . 

Even  with   a  large   computer  system,    individual  working 
space   is  not   usually   large  enough   for   storing  all   data 
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Pixel  Value 


Object-Coordinate 


Viewing -Coordinate 


Ob j  ec t -Coordinate 


(A) 


(B) 


Figure  3.10 


Interpolation  Kernel  over  ob iect -  coordinate (A) 
and  Bilinear  Int erpolat ion(B ) . 


values.    Therefore   this   method   is   not     practical   for 
manipulation  of  3D  data. 

A  displayed  image  obtained  by  this  method  is 
shown  in  Figure  3.11.  It  shows  a  better  result  than  the 
interpolation  after  rotation  of  obj ect -  coordinate .  Blurring 
effect  is  due  to  inherent  linear  interpolation  error. 

c.   Interpolation  Using  Cone-Shape  Kernel 

We  have  already  disscussed  two  interpolation 
methods.  But  these  two  methods  have  flaws  which  make  it 
difficult  to  apply  it  to  the  3D  situations.  The  third 
method  which  avoids  these  flaws  uses  a  2D  cone-shape  kernel. 
Figure  3.12  illustrates  how  the  cone-shape  kernel  is  used 
for  2D  linear  interpolation.  The  density  distribution  of  a 
pixel  in  the  exact  2D  linear  interpolation  should  be  done 
like  Figure  3.12(B).  But  the  cone-shape  density  distribution 
method  is  used  for  this  purpose.  The  advantage  of  using  a 
cone-shape  kernel  is 
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(A)  Original  Image 


(B)  15  Degrees  Rotation 


(C)  30  Degrees  Rotation 


(D)  45  Degrees  Rotation 


Figure  3.11    Displayed  Image  by  Interpolation  of  the  Signal 
from  the  Nearest  4  Samples  at  the  Grids 
of  the  Ob j ect -  Coordinate . 


(1)  It  is  invariant  to  the  direction  of  rotation. 

(2)  Interpolation  is  performed  on  the  obj ect -  coordinate 
The  algorithm  for  this  method  is 

(1)  Consider  one  obj ect -  coordinate  grid. 


46 


Pixel  Value 


Pixel  Value 


Object -Coordinate 

(A)  (B) 

Figure  3.12    (A)  Cone-Shape  Kernel  (B)  Rectangular  kernel. 

(2)  Transform  it  to  the  viewing- coordinate  value  by 
coordinate  rotation. 

(3)  Calculate  the  distance  from  it  to  the  4  surrounding 
obj ect -  coordinate  grids. 

(4)  If  a  distance  is  smaller  than  unity,  calculate  the 
distribution  component  of  a  intensity  value  inversely 
proportional  to  the  distance.  Otherwise  assign  zero. 

(5)  Sum  the  distributions  of  all  points  to  get  the 
intensity  value  at  the  grid  of  the  object  coordinate. 

Using  this  kernel,  only  viewing- coordinate  grids  inside 
the  kernel  are  affected  by  the  sample  value  at  the  center  of 
the  cone.  Otherwise  it  is  ignored.  Figure  3.13  indicates 
that  the  error  introduced  by  using  the  cone-shape  kernel  is 
not  severe.  Resulting  displayed  image  is  almost  the  same  as 
the  above  one.  Because  this  cone-shape  kernel  is  invariant 
to  the  direction  of  rotation,  it  can  be  extended  easily 
into  a  sphere-shape  kernel  in  a  3D  situation. 
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(A)  Original  Image 


(B)  15  Degrees  Rotation 


C)  30  Degrees  Rotation 


(D)  45  Degrees  Rotation 


Figure  3.13    Displayed  Images  by  Interpolation 
Using  Cone-Shape  Kernel. 


The  disadvantages  of  this  method  are:  (1)  It 
involves  many  multiplications  to  calculate  the  distances. 
(2)  Error  includes  inherent  linear  interpolation  error  plus 
the  cone-shape  kernel  error.  The  advantage  of  this  method 
is  that  it  can  easily  be  applied  in  a  3 -dimensional 
situat  ion . 
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3 .   Three -Dimensional  Interpolation 

As  dicussed  in  the  2D  situation,  linear  interpola- 
tion using  a  cone-shape  kernel  is  the  most  appropriate 
method.  Because  the  shape  of  this  kernel  is  invariant  with 
respect  to  the  direction,  it  can  be  easily  extended  to  the 
3D  interpolation.  Applying  the  linear  interpolation  method 
to  the  3D  situation  requires  a  very  complicate  procedure. 
After  rotation,  a  parallelepiped  looks  like  a  polygon 
instead  of  a  cube  to  the  sight  of  the  viewer.  Therefore, 
performing  linear  interpolation  to  the  sphere  is  easier  than 
to  the  polygon.  In  this  case,  only  grids  of  object- 
coordinate  inside  of  the  sphere  are  affected  by  a  pixel 
value  lying  at  the  center  of  a  sphere.  The  algorithm  for 
this  method  is  : 

(1)  Consider  one  obj ect -  coordinate  grid. 

(2)  Transform  it  into  viewing- coordinate  value  by 
coordinate  rotation. 

(3)  Calculate  distance  from  it  to  the  surrounding  8 
obj ect -  coordinate  grids. 

(4)  Calculate  the  distribution  of  intensity  value  at  the 
8  obj ect -  coordinate  grids  when  a  distance  is  smaller  than 
unity.   Otherwise  ignore  it. 

(5)  Sum  the  distribution  components  of  the  8  surrounding 
intensity  values  to  calculate  the  value  at  a  object- 
coordinate  grid. 

The  disadvantages  of  this  method  are  (1)  Necessity  for 
many  floating-point  multiplications,  (2)  Total  error  is  the 
sum  of  inherent  linear  interpolation  error  and  sphere-shape 
kernel  error.  More  details  will  be  discussed  in  the  Chapter 
IV. 
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IV.  PROGRAM  IMPLEMENTATION  OF  REPROJECTION  AND  DISSECTION 

A.   ORTHOGRAPHIC  REPROJECTION 

Orthographic  repro j ection  is  performed  numerically  in  the 
computer  by  summing  the  value  of  voxels  along  parallel  paths 
through  the  reconstructed  volume.  The  reprojected  image 
constitutes  a  2-dimensional  image  on  the  screen. 

Figure4-1   illustrates   the    orthographic   reprojection 
process  onto  the  X*-Y#  plane  along  Z»  coordinate. 


Reprojected 
Pixels  on  the 


(Yk  th  Row 


Viewing-Plane  ^       • L__S»^Lr4 


Viewing-Plane 


Summation  Paths 
(Yk  th  Slice  of 
Volume ) 


Figure    4.1         Reprojection    Process    before    Rotation. 

To       implement    this      process       in      the    program,  first       a 

viewing-plane    (2-dimensional    array:  'Viewpln')       is    defined 

parallel  to  the  viewing- coordinate  X»-Y»  shown  as  in  program 
'Rot_3_Dim' .  The  viewing-plane  composes  of  X* (column)  and 
Y* (row)  coordinates  and  lies  orthogonally  to  Z •- coordinate 
(X#,Y»,Z#    are    viewing- coordinates  )  .  The    Yic'th    slice    of    the 

volume  is  projected  onto  the  Yk  ' th  row  of  pixels  in  the 
viewing-plane    as    shown    in    Figure    4.1. 
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Before  rotation,  every  voxel  is  arranged  along  viewing- 
coordinate  so  obj ect -  coordinate  grids  correspond  to 
viewing- coordinate  grids. 

The  reprojection  can  be  performed  as  in  the  following 
procedure : 

(1)  For  every  incoming  voxel,  read  the  viewing- 
coordinate  Xft» ,  Yft* ,  then  sum  the  density  value  onto  a  pixel 
corresponding  at  X^#,Y|^#  in  the  viewing-plane  (  array: 
'Viewpln' ) . 

Viewpln(XK« ,Y^-)  :=  Viewpln(XK • , YK* )  +  Pixel (XK • , YK • , ZK* ) \ 

(2)  Continue  to  reach  the  last  voxel. 

After  rotation,  it  is  not  clear  what  portion  of  a  voxel 
will  contribute  the  pixel  value  in  the  viewing-plane  because 
each  voxel  looks  like  a  polygon  to  the  sight  of  view. 
Through  the  interpolation  method  discussed  in  the  Chapter 
III,  the  appropriate  portions  of  the  voxel  value  are  added 
to  the  pixel  values  in  the  viewing-plane  (see  subroutine 
'Projection'  and  'Dissection'  of  the  program  'Rot  3  Dim'). 


Re projected 
Pixels  on  th 
Viewing-Plan 


Viewing-Plane 


Summation  Paths 


Figure  4.2    Reprojection  Process  after  Rotation 
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The  resulting  pixel  values  in  the  viewing-plane  will 
exceed  the  allowed  pixel  values  of  the  display  device  so 
that  the  pixel  should  be  rescaled  to  fit  the  specific 
display  device  requirement . 

Several  reprojected  images  with  sphere-shape  kernel  are 
illustrated  in  the  Figure  4.3.  Different  viewing- angles 
provide  more  informations  about  the  structures  of 
3-dimensional  data. 
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30  Degrees  45  Degrees 

(A)  Rotation  about  X-Axis 


75  Degrees 


90  Degrees 
(B)  Rotation  about  Y-Axis 


180  Degrees 


30-30  Degrees        45-45  Degrees        60-60  Degrees 
(C)  Rotation  about  Both-Axis 

Figure  4.3    Displayed  Images  by  Reproj ect ion . 
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B 


SINGLE  PLANE  DISSECTION 


The  dissection  is  the  technique  that  an  object  plane  is 
clearly  displayed  showing  spatial  relationship  with  other 
planes  in  the  volume  data.  But  only  single  plane  dissection 
is  considered  in  this  work  because  the  dissection  process  is 
very  similar  to  the  single  plane  dissection.  Single  plane 
dissection  is  performed  by  removing  all  other  planes  and 
displaying  an  object  plane. 

Suppose  we  want  to  see  a   plane  to  the  selected  viewing- 
angle  after  rotation.   An  viewing-plane  will  be  lying  in  the 


'  . Viewing-Plane 


Figure  4.4    Single  Plane  Dissection. 

Viewing- coordinate  shown  as  in  Figure  4.4.  After  rotation, 
viewing-plane  cuts  the  volume  with  an  arbitrary  angle.  The 
voxels  reprojected  onto  the  viewing-plane  are  not  squares 
rather  polygons  with  an  arbitrary  shape.  Therefore  it  is 
difficult  to  determine  what  portion  of  sample  values  will  be 
contributed  to  the  pixels  in  the  viewing-plane. 

To  calculate  the  intensity  value  of  pixels  in  the 
viewing-plane,  3 -dimansional  interpolation  is  necessary.  The 
interpolation  using  sphere-shape  kernel  is  used  to  calculate 
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the  intensity  values  on  the  viewing- coordinate  grids. 
Subroutine  'Dissection"  in  the  program  'Rot_3_Dim'  shows  the 
dissection  process.   Procedure  for  this  scheme  is 

(1)  Accept  a  desired  rotation  angle  and  desired 
Z  •- coordinate  to  see. 

(2)  Take  a  voxel  from  the  data. 

(3)  Transform  it  into  the  corresponding  viewing- 
coordinate  value. 

(4)  Devide  the  voxel  value  to  the  surrounding  8 
coordinate-grids  by  3-dimensional  interpolation  method. 

(5)  If  the  Z» -  coordinate  of  the  grid  corresponds  to  the 
desired  Z •- coordinate ,  add  the  portion  of  the  voxel  value  to 
the  corresponding  pixel  in  the  viewing-plane .  Otherwise 
ignore  it . 

(6)  Continue  to  the  last  voxel. 

Several  dissected  planes  from  the  reprojected  image  are 
illustrated  in  Figure  4.5.  We  can  see  the  checker-board 
effect  in  some  displayed  images  which  does  not  appear  in  the 
figure.  This  means  that  the  linear  interpolation  with 
sphere-shape  kernel  is  not  sufficient  for  the  single  plane 
dissection  in  the  certain  rotation  angle.  If  one  want 
higher  quality  images,  more  precise  and  efficient 
interpolation  method  should  be  developed. 
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45  (ZCRD=25)  60  (ZCRD--25)       90°(ZCRD=0) 

(A)  Rotation  about  A-Axis 


45  (ZCRD=20) 


135W(ZCRD=20)        80°(ZCRD=0) 
(B)  Rotation  about  Y-Axis 


30  -6CT  (ZCRD=0)  45  -45  (ZCRD=0) 

(C)  Rotation  about  Both-Axis 

Figure  4.5    Displayed  Image  by  Dissection 
(ZCRD  is  the  Z  th  Slice  of  Volume  Image). 
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V.  CONCLUSION  AND  RECOMMENDATIONS 


A.   CONCLUSION 

An  inherent  disadvantage  of  the  reprojection  method  is  the 
obscuring  of  underlying  and  overlying  structures  in  the 
2-dimensional  viewing- screen.  Therefore  it  is  desirable  to 
have  better  visibility . of  selected  portions  of  the  volume 
data.  Interactive  user  selection  of  an  adequate  viewing- 
angle  and  the  application  of  the  dissection  method  is 
proposed  to  resolve  this  problem.  Another  important 
problem  is  the  need  to  manipulate  a  huge  amount  of  data 
as  efficiently  as  possible. 

Three  different  algorithms  were  studied  here  for  the 
coordinate  transformation:  Direct  matrix  multiplication, 
Decomposition  with  fixed  coordinate,  and  Recalculation  and 
Indexing . 

In  the  Precalculation  and  Indexing  method,  the  compo- 
nents of  coordinate  and  matrix  element  are  calculated  itera- 
tively  and  stored  in  an  array.  Then  these  components  are 
indexed  according  to  the  incoming  coordinates  and*  used  in 
the  calculation  so  that  the  number  of  matrix  multiplication 
is  decreased.  This  method  requires  an  appropriate  size  of 
memory.  It  is  possible  to  rotate  643  data  elements  in  25 
seconds  with  the  VAX-11/750  minicomputer. 

To  realize  this  algorithm  in  hardware,  two  hardware 
design  schemes  have  been  suggested  in  Chapter  II. 

After  rotation  to  a  selected  direction,  the  positions  of 
the  voxels  do  not  correspond  to  the  viewing- coordinate 
grids.  Therefore  several  linear  interpolation  algorithms  are 
studied . 

The  reason  why  linear  interpolation  using  a  cone-shape 
kernel  can  be  extended  from  2D  to  3D  is  attributed  to  its 
invariant  shape   with  respect   to  the   different  directions. 
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The  sphere-shape  kernel  is  a   3D  extension  of  the  cone-shape 
kernel  in  the  2D  situations.  p 

Using  this  3D  interpolation  method,  a  reprojection  and 
single  plane  dissection  program  has  been  tested  with  an 
artificial  data.  Although  linear  interpolation  using 
sphere-shape  kernel  has  inherent  errors,  this  method 
provides  the  reasonable  quality  images  in  a  reasonable 
amount  of  time. 

Manipulation  of  a  huge  amount  of  data  such  as  3D 
computed  tomography  data  requires  a  continuous  trade-off 
between  processing  time  and  memory,  or  between  processing 
time  and  image  quality.  The  only  way  to  resolve  this  kind 
of  problems  is  to  use  heuristic  approaches.  First,  select 
an  initital  algorithm  and  test  it.  If  the  result  is  not 
adequate  for  the  specific  purpose,  try  another  algorithm. 
This  procedure  is  continued  until  an  appropriate  result  can 
be  reached. 

B.   RECOMMENDATIONS 

The  coordinate  rotation  and  reprojection  method  can  be 
widely  used  in  the  viewing  or  manipulating  of  3 -dimensional 
data.  The  algorithms  which  have  been  developed  in  this  work 
may  not  be  optimum  because  these  are  the  result  of  heuristic 
approaches . 

The  coordinate  transformation  algorithm  implemented  here 
can  rotate  any  reasonable  size  of  3D  data  in  a  short  time. 
But,  further  studies  are  required  in  order  to  get  better 
quality  images  from  a  3D  interpolation  algorithm. 

It  is  worth  while  to  test  the  implementation  of  the 
reprojection  algorithm  with  real  CT  scan  data.  In  this  case 
new  problems  may  arise,  which  will  require  further  studies. 

The  implementation  of  the  dissection  and  the  dissolution 
capabilites  are  very  helpful  for  the  investigation  of 
structures  in  3D  volume  data. 
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